


(HASA-CR-3169) A SIMPLIFIED COHPOTER 
PROGRAM FOB THE PREDICTION OF THE LINEAR 
STABILITY BEIIAVIOB OF LIQOID PROPELLANT 


COHBOSTOBS (Colorado State Univ.) 39 p 

HC A04/BF A01 CSCL 21H H1/20 


N79-28226 

Onclas 

34097 








NASA Contractor Report 3169 


A Simplified Computer Program 
for the Prediction of the 
Linear Stability Behavior of 
Liquid Propellant Combustors 


C. E. Mitchell and K. Eckert 
Colorado State University 
Fort Collins, Colorado 


Prepared for 

Lewis Research Center 

under Grant NGR-06-002-095 


fW\SA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Branch 


1979 



TABLE OF CONTENTS 


NOMENCLATURE v 

INTRODUCTION 1 

THEORY 3 

Combustor Model 3 

Method of Solution 10 

COMPUTATIONAL METHODS 15 

Matrix Sizing and Program Convergence . . 15 

Program Options 18 

Choice of Fundamental Acoustic Mode ... 21 

REFERENCES 23 

Program MODULE 25 

General Description of Program .... 25 

Program Input . 32 

Program Output 34 

Sample Run 34 

Program MODULE Listing 40 


iii 



NOMENCLATURE 


Letters 

a - speed of sound 

AMF - aperture mean flow 

BR - backing distance for slot absorber 

F - quantity for integral ^verning equation 

f - quantity for integral governing equation 

6^ - modified Green's function, defined after Equation 15 

i - / -l 

J - Bessel Function of first kind of order m 

m 

K - acoustic impedance of slot absorber 

L - nondimensional chamber length 

i - radial acoustic mode assumed 

L, - actual length of aperture for slot absorber 

a 

L g ^ - effective length of aperture for slot absorber 

M - mean flow Mach number 

m - transverse acoustic mode assumed 

m - mass generation rate 

n - longitudinal acoustic mode assumed 

n - interaction index 

n - outward directed normal unit vector 

P - pressure 

r - radial dimension 

r - radius of chamber 

c 

Ro - resistance of slot absorber 

S - surface of combustion chamber over which integration is 

to be carried out 
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T - temperature 

t - time 

u 0 - normal component of velocity osci.lation in 

c tangential direction 

u - normal component of velocity oscillation in 

r radial direction 

u t - transverse veloc. -y 

V - velocity or volume of combustion chamber 

Wa - aperture width for slot absorber or length 

of acoustic liner 

X - distance from injector face to beginning of 

a slot absorber or acoustic liner 

X, - distance from injector face to end of slot 

° absorber or acoustic liner 

z - longitudinal dimensi . 

Greek Letters 

6 - acoustic admittance of a surface 

y - ratio of specific heats 

e - nondimensional wave amplitude 

n - acoustic eigenvalue 

e - angle in radians 

A - normalization factor defined after Equation 15 

X. - root of Bessel Function of first kind, such that 

*" j ;<v ■ 0 

y - coefficient matrix 

p - density 

t - sensitive time lag 

$ - velocity potential 
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H< - normalization factor defined after Equation 15 
" acoust ’ c eigenfunction 
u) - complex frequency 

Superscripts 

-*■ - vector quantities 

* - dimensional quantity 

' - derivative with respect to argument, or perturbation quantity 

- mean or steady state quantity 

Subscripts 
I - injector 

L - liner 

N - nozzle 
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INTRODUCTION 


The purpose of this report is to present an analytical technique and 
a computer program which can be used for the prediction of the linear 
stability behavior of liquid propellant combustors. The technique involved 
has been developed over the last few years at Colorado State University ir> 
the examination of several aspects of the instability problem. Basically, 
the approach employs a Green's function integral method in the iterative 
determination of combustor frequency, decay rate and spatial waveform. 

This general approach has been applied to several different combustor 
models in the examination of different aspects of the linear instability 
problem. (Ref. 1-10.) This work was performed by several different people 
(mainly graduate students), and a wide variety of nomenclature and program- 
ming techniques has resulted. The details of the analytical approach have 
also varied from author to author though the general method remained the 
same. This more or less comprehensive compendium of programs and analyses 
as it exists in its several forms is cumbersome, somewhat redundant, and 
certainly hard to use as a designer's tool. 

With this in mind it was decided to develop a simplified stability 
analysis and computer program which contained the most important features of 
the earlier work in a format that would be relatively easy to use. Conse- 
quently, the main goal of this effort has been the development of a computer 
program simple enough to be used effectively by a person without an exhaustive 
background in either advanced mathematics or stability theory. 

In order to do this some compromises have had to be made as far as 
comprehensiveness and accuracy are concerned, and some aspects of the stability 


problem treated previously have not been included. For example, the effect 
of distributing combustion sources along the combustor axis (as opposed to 
having a concentrated combustion zone near the injector) on overall stability 
has been studied and analyzed using two different approaches (Ref. 8, 9). 

This effect is not included in the simplified model presented here, however. 

I 

The justification for this is based on the fact that much greater complexity 
is introduced into both the analysis and the computer program when distrib- 
uted sources of combustion are considered, while the qualitative stability 
behavior is very similar to that predicted for concentrated combustion. 
Moreover, the quantitative effect of distributing the combustion is stabi- 
lizing relative to the predictions for a concentrated combustion zone. Thus, 
the simplified model presented here will tend to give conservative estimates 
of combustor stability when the combustor being examined has its combustion 
zone well distributed (axially). 

Other effects such as irrotational ity, entropy variations, and droplet 
drag effects have also been ignored since their influence has been found to 
be small, stabilizing or both. 

The body of the report is divided into three main sections. The first 
(called "Theory") presents the model and method or analysis. The second 
section (called "Computational Methods") presents the basics of the computa- 
tional method and the user options available. The final section (called 
"Program MODULE") gives a user's manual, sample input and output and a flow 
chart. It is not necessary for a person wishing to use the computer program 
(MODULE), to follow the analytical details of the first section. It will be 
necessary, however, for him to understand the basics of the model and general 
method of approach as presented in that section so that appropriate input to 
the program may be made and correct interpretation of the output can result. 
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THEORY 


Combustor Mode l 

The motor configuration considered here Is characterized by circular 
cylindrical geometry, a concentrated combustion zone located at the injector 
end of the combustor, a nozzle at the opposite end, and either an acoustic 
liner or a slot absorber located in the cylindrical walls. A sketch of 
the combustor model is given in Figure 1. In the development of a linear 
stability model for a combustor of this type it is first necessary to repre- 
sent the four main features of the configuration using appropriate mathemat- 
ical models. The four aspects of the problem requiring such modeling are 

1) The gasdynamic flow field 

2) The combustion zone 

3) The nozzle 

4) The acoustic liner or absorber. 

Each of these will be discussed separately before going on to a presen- 
tation of the global stability model and analytical technique. 

1) The gasdynamic flow field 

The flow field downstream of the concentrated combustion zone is taken 
to consist of a single component, single phase product gas which is non- 
conducting, inviscid and calorically perfect. The flow is assumed to be 
homentroplc and irrotatlonal . As long as the combustion zone is concentrated 
and pressure waves are of small amplitude, it has been shown that these 
approximations are not severely limiting and self-consistent (Ref. 11, 12). 
Before presenting the equations describing this flow field the relevant 
state and flow variables are non-dimensional 1 zed as follows. 
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P 3 p*/p* , T 3 T*/T* 

V 3 V*/a* . P * P*/P* 

The independent variables are nondiinenslonalized as follows. 

t = t*(r */a*) , r = r*/r* 

w L 

z = z*/r* 

where * denotes dimensional quantities and ~ denotes mean chamber values. 
Using this nondimensional scheme the cr servation equations become 


22. 

dt 

+ 

V»(pV) = 0 

CONTINUITY 

(1) 

DV 

Pot 

+ 

-* 

- VP = 0 
Y 

MOMENTUM 

(2) 

p * 

P Y 


HOMENTROPIC 

(3) 

p = 

pT 


STATE 

(4) 

V * 

-*■ 

v$ 


IRROTATIONALITY 

(5) 


Under the assumption of small amplitude oscillations the state and flow 
variables are represented as the sum of a mean (steady state) component and 
an oscillatory component, products of which are ignored as being higher order 
terms. Thus 

P=l+p' $ = Mz + <J> 

P = 1 + P 

T = 1 + T' V * Me z + V<J> 
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where e z is the unit vector in the axial direction and M is the mean flow 
Mach number. 

After some manipulation the conservation equations can be reduced to a 
simple scalar partial differential equation 

f2 0 - 0 * 2M lift + * U (6) 


The equation relating the state variables to 4 > is 
p... Y( |i*„||) (7) 


Periodic oscillations in time are assumed so that 

p' = p(r, ez) e 1u)t 
4> = <D(r, 0, z) e iwt 

where to = to R + iX is the complex frequency, to^ the frequency, 

X the decay rate. (If X > 0 decay occurs.) 

2) Combustion zone response model 

It is assumed that all combustion occurs in a length small compared with 
the combustor's axial dimension. In the steady state mass is produced at the 
rate m = M in the nondimensional system used here. No attempt to describe 
the details of the combustion process is made. Instead it is simply assumed 
that the combustion zone is sensitive to pressure oscillations and responds 
to these oscillations through a combustion zone admittance function Bj. 

Thus 


V<J> • n g p' (z = 0) 


( 8 ) 
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f 


Bj is taken to be a constant for t’e entire combustion zone, though p' 
is, of course, a function of r and z as well as time- In terms of the 
mass perturbation rate, m' , the response condition is 

- (” - 8,) P' (9) 

Sj is, in general, complex so that all phasings between m' (or u') and 
p' are possible. Note that if the real part of $j is greater than — , 
the combustion zone provides a damping rather than driving effect. 

It is also possible to relate Bj to the interaction index n, and 
time lag t of the Crocco sensitive time lag model. The appropriate 
relationship is 

8, ■ - n (1 - e _i “ T ) ) . (10) 

Values of Bj (or n, and t) must be supplied by the program user or 
calculated as output, given all other parameters. These options will be 
discussed later. 

3) Nozzle model 

Here again no attempt is made to investigate the details of the nozzle 
flow and, instead, a nozzle admittance function B N is used. 

V<f> • n = B f , p‘ (z = L) (il) 

Values of 6^ are to be supplied by the user. Tables of admittance 
functions are given in Ref. (13), for example, for conical nozzles. In the 
absence of any knowledge of the nozzle response value it is suggested that 
the simple "short" nozzle value 



be used. 
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4) Acoustic liner or slot absorber model. 

Two possibilities are considered. The first is an acoustic liner of 
uniform average admittance, 6^ , which is uniform in the azimuthal (6) direc- 
tion and extends along the cylindrical wall from z = to z = Xg . For 
this liner the appropriate boundary condition is 

-V 

V<J> • n = $ L p' (12) 

No attempt is made to calculate 8^ in either the analysis or computer 
program and therefore B L 'xist be supplied by the user. 

The second absorber configuration considered is a circumferential slot 
machined into the cylindrical wall of the chamber and acting as a Helmholtz 
resonator. The geometry assumed is shown in Figure 2. All dimensions are 
nondimensional through division with the chamber radius. 

The appropriate boundary condition at r = 1 (chamber wall) over the 
aperture width W A (=[x B -x A ]) 

V<j> • n = 8 l P 

or $j> • n = p 

where K = - 4 - is the impedance at the aperture entrance. K is used in 
YP l. 

this case to be consistent with existing treatments of Helmholtz resonators 
of this general type. 8. and K are, in general, complex with K = R + ik , 
where R Q is the resistance, k the reactance. Standard relationships for 
R q , resistance and taken from Reference (15) and (16) respectively, 
are given below. 
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Figure 2 . Slot Absorber Geometry 
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. 8^1 . 


5 ( AMF ) R q + 



(15) 


where e is wa ve amplitude and p is the modulus of the pressure at r = 1. 


L eff = L A + (0-375)(0.85 )W a 


1 - 0.71 


w aV s 


(5) 


An expression for the reactance k, comes directly from linear Helmholtz 
resonator theory and is given by 


k ~ W R p ap L eff 


W A ^a p a 
to R W^. BR 


where p 

ap 




are, respectively. 


the nondimensional aperture density, absorber cavity mean sound speed, and 
absorber cavity density, and oj r is the real part of the oscillation frequency. 
All these quantities, as well as the geometrical quantities in Fig. 2, AMF, 
and the assumed wave amplitude, c , must be supplied by the user. The 
expression for R q , Equation (4), is then solved iteratively for R Q as a 
function of frequency. Thus, an expression for K(uj) (or B l (oj)) is found 
numerical ly. 


Method of Solution 

The governing partial differential equation. Equation (6) along with the 
necessary boundary conditions (Equations (8) (10) (12) or (13)) are transformed 
to integral form using a Green's function, and the resulting integral equations 
are solved iteratively. .Details of the transformation and solution method are 
presented in References (1, 2, 7). Only those relationships, definitions, and 
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equations necessary for understanding and using the computer program which 
determines combustor stability will be presented here. 

The transformed integral equations for 4>{r, 0, z) and 

u) are 


4> 


"f 

£mn 


///S N (r/r o ) F ] (4))dV o 


+ // G N (?/r D ) f 1 (4>)dS o 

S o 


(14) 


0) 


- nL, = ///ft — F, (4>)dV + //fi A _ f , ( ij) ) dS 
S.mn „ imn 1 ' im 1 


(15) 


where F, (4>) = 2icoM f|- + M 2 

• oZ dZ 


fi » - 0P 


3 = 3^ at nozzle (z = L) 

3 = Bj at combustion zone (z = 0) 

3 = \ (or ^ ) at liner (or absorber) (r 
3 a 0 on all other surfaces 

G N (r/r Q ) = ZZL — m - 7 - ( y - 


= 1 ) 


Jimn 


' "Im > 


S. f £ . m j* m , n^n, simultaneously 



- 12 - 


J(L m r) cos — - cos m 6 
n m Jem L 

£mn A 

A £mn 

A lmn ' W [ J m (X Jtm r) c05 cos m e ] 2 dV 

x <tm are the roots of J' m (A^ ) = 0 



^£mn are the normalized eigenfunctions for a cylindrical chamber with 
no mean flow and non reactive walls. I , m, n are the set of integers 
giving the radial, azimuthal, and axial character of the particular eigen- 
function (or acoustic mode) in question. Thus, represents a first 
transverse mode, fi-^g a secor 'd transverse mode, J^qo a first radial mode, 
ftggi a first axial mode, a combined first transverse first axial mode, 
etc. The associated eigenvalues (acoustic frequencies) are 



The solution technique revolves around the assumption that the actual 
solution including mean flow and reactive walls has a character that is 
reasonably close to one cf these acoustic modes. The particular acoustic 
mode most characteristic of the overall oscillation is called . where 


tmn 

i , m, n are the associated indices giving the radial, azimuthal, and axial 
character. The related eigenvalue (acoustic frequency) is njftn • A discussion 
of the selection of in applications will be given later. 
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The equation for <j> , Equation (14)^ implies that <p takes the following 

form 


* in 

£nm 


limn 




cos 


nrrz 


cos m 6 


where the coefficient matrix u^ mn is determined by evaluation of the 
integrals on the right hand side of Equation (14). 

Because of the symmetry in the 0 direction which results from the 
assumptions concerning the boundary conditions, the series in m actually 
contains only one term, m . 

Thus, for the model used here $ may be written 






(16) 


2 IT 

where e 2 = /(cos m e) 2 d0 . Exactly the same coefficient matrix would 
9 o 

result if traveling waveforms were assumed. In this case 

* = *(r, 2 )e 1( " t + * 0 > 

* • f % J s <V> “> s “r (,7) 

and would be identical to the standing wave matrix. 

The matrix p £n and the complex frequency to are determined by an 
iterative process. The lowest order guess for <J> (or u^ n ) is used in 
the integral expressions of Equations (14) and (15) to compute improved 
values for the p^ n and to . The process continues until successive 
iterations are invariant to some degree of accuracy. A natural choice for 
the lowest order estimate for $ would be ; the lowest order frequency 

would then be • Though these initial guesses will work in general, 
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experience has indicated that convergence can be slow and matrix sizes 
large, particularly when the mean flow Mach number is greater than about 
0.3. Better convergence and a smaller matrix size are possible if the 
separation of variables solution for a combustor with mean flow but without 
an absorber is used. This solution was originally developed by Priem and 
Rice (Ref. (14)); in the modified form appropriate here it is discussed in 
References (1) and (2). The computer program presented later uses this form 
as the lowest order <f> . 

In addition to assuming a lowest order form for <p and oj and iterating, 
it is also possible to fix co at some prescribed value (supplied by the user) 
and iterate to find the appropriate and B T (or n, and r ) from the 

same equations. The latter approach is used to solve for the combustion 
response necessary to sustain an oscillation of a given frequency and decay 
(growth) rate and known absorber and nozzle admittances. It would also be 
possible to set up the technique to solve iteratively for another parameter, 
such as nozzle admittance, for given combustion admittance; however, this has 
not been done in the program presented here. 
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COMPUTATIONAL METHODS 


As discussed in the "Theory" section. Equations (14) and (15) are 
set up for iterative solution. Computer program MODULE is an algorithm 
for performing the necessary iterative computations on a digital computer. 
Several different choices are possible as far as input, output, and 
accuracy are concerned. These choices will be discussed in this section. 


Matrix Sizing and Program Convergence 

In the solution of Equations (14) and (15) two variables are always 
iterated. One of these is the perturbation velocity potential <J >(r,0,z). 
The other is either the complex frequency oj (co^ + iA) or the complex 
combustion admittance, Bj (real ( Bj ) + i imag (Bj)). 

The perturbation velocity potential is represented by a series 
expansion (Equation (16)). 


4>(r,0,z) 


II 

£n 


£n m 


<\m r) 


COS 


nuz 

L 


cos m 8 

e 0 


Thus, solution for the coefficient matrix yields <j>( r, 0, z ) and, 

in fact, it is this matrix which is the actual iterated variable in the 
solution algorithm. Formally, is doubly infinite in £ and n . 

That is, lS£<°°,oSn<«. As a practical matter, however, limits 
on the largest values £ and n may take (in other words the dimensions 
of matrix M^ n ) must be determined. It should be recalled here that the 
integers n are associated with axial dependence (through cos ^p) while 
the integers £ are associated with radial dependence (through (A^r)). 

Any choice for the maximum number of "£ terms" and "n terms" will 
limit accuracy. A compromise between program run time, storage requirements. 
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and accuracy is desirable. Naturally, no one choice will be optimal for 
all combustor configurations. However, hundreds of runs with "typical" 
designs have indicated some rules of thumb to be used. 

First of all, in none of the combustors investigated was any signif- 
icant increase in accuracy obtained by keeping more than 50 terms in the 
axial direction or ten terms in the radial direction. That is, keeping 
100 terms in the axial direction or 20 terms in the radial direction 
affected the values of the iterated variables only very slightly (<0.25%). 
Consequently, the program as written accepts a 10 x 50 matrix size for 
as the maximum allowable. In the program variables this means LTS - 10, 

NTS ^ 50, where LTS and NTS are, respectively, the number of terms in the 
direction and the number of terms in the "n" direction. 

The question as to the "best" values of LTS and NTS to use in a given 
combustor configuration is difficult to answer. Eckert (Ref. (14)) has 
studied optimal values for LTS and NTS for a "typical" configuration and 
suggests values of 3 for LTS and 16 for NTS. However, for a combustor with 
no absorber, a single term (l) is necessary for description of the radial 
field and LTS = 1 in this case. On the other hand, if the Mach number is 
small, mean flow effects are less important and fewer terms in the axial 
direction (smaller NTS) would be needed. However, for configurations with 
large absorber effects or high Mach numbers (> 0.4) it is likely that "best" 
values for LTS and NTS could be greater than 3 and 16, respectively. 

With this in mind it is suggested that the values LTS = 3 and NTS = 16 
be used as a general rule. If strong absorber or high Mach number effects 
are present and may compromise accuracy, it is suggested that results with 
LTS = 9 and NTS = 50 be computed and compared with the smaller matrix 
results to estimate accuracy. Values of LTS and NTS larger than 3 and 16 
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could then be inserted until the desired accuracy relative to the 9 x 50 
size was obtained. It should be noted that improvement in accuracy is 
monotonic with increasing NTS. The same is not true for LTS because of the 
alternating nature of the series involved and best results occur if LTS is 
an odd number (3, 5, 7, 9). 

Once the dimensions of the p^ n matrix are determined it is next 
necessary to decide upon an acceptable convergence condition for the itera- 
tion process. The second iterated variable (either ou or gj depending 
upon the application) is used to do this. Successive values of the iterated 
variable are compared. When the difference between the two values is less 
than some value, adequate convergence is assumed. Since both ^ and gj 
are complex numbers, it is necessary that both the real and imaginary parts 
converge in the sense just mentioned. In this program, however, it is 
convenient instead to deal with w (or gj) in complex polar notation, and 
require that successive values of the modulus and phase angle converge. 

This is because the phase angle is frequently rear zero and can cause 
problems in the definition of convergence for the imaginary part of the 
iterated variable. In program MODULE convergence is assumed when the percent 
change in the modulus of the iterated variable is less than the value ERROR 

and, at the same time, the absolute value of the change in the phase angle 

-5 

is also less than ERROR. ERROR can take values between 10 and 1.0. 

In most cases convergence to within 0.1% or less is rapid, usually 
occurring in ten iterations or less. However, for some choices of para- 
meters and for some program options it can be much slower or not occur at 
all. For this reason a maximum desired number of iterations must be spec- 
ified. This is done through program variable IDMAX which can take any 
integer value. If convergence does not occur in the number of iterations 
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specified by IDMAX, the iterative loop terminates, and program values at 
the last iteration are output. 

Program Options 

In addition to choosing either cj or Bj as the iterated variable, 
choices are possible as far as the form of the combustion response model 
and the acoustic absorber. Taken togetner this results in six distinct 
ways of running the program. These possibilities are labelled options and 
are described sequentially below. For all of the options it is necessary 
that certain design or program variables be specified by the user. These 
parameters are y (ratio of specific heats), M (mean flow Mach number), L 
(chamber length to radius ratio), 8^ (complex nozzle admittance), ERROR 
(maximum error allowable in determining convergence), ITS (number of terms 
in radial direction kept), and NTS (number of terms kept in the axial 
direction). 

Option 1 

This option is designed to compute frequency and decay rate (complex 
frequency) for known combustion zone admittance and known acoustic absorber 
(or liner) length and admittance. The iterated variable is the complex 
frequency. Required to be input to the program are 8j , B L , X A and X g . 

X^ and X g are the nondimensional distances to the start and end of the 

acoustic absorber, respectively. Output are co R and X, , the input 
parameters, and n and t , the interaction index and time lag correspon- 
ding to the given and the converged value for oj . 
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Option 2 

Option 2 is similar to Option 1 except that the <ombustion response 
is described by n and x instead of 3j . In this case n and x are 
input and Bj is output. Other input and output parameters are the same 
as for Option 1 . 

Option 3 

In this option the acoustic absorber is of the slot design type des- 
cribed earlier. The combustion response is described by 3j , and u> is 
the iterated variable. Required input variables are Bj , B^ (absorber 
backing distance), kl (absorber cavity width), L, (absorber aperture length), 
X A ’ X B ’ ®a ’ (ratio of sound speed in the cavity to sound speed in the 
main chamber), p, (nondimensional aperture gas density), AMF (aperture mean 

a 

flow), and e (wave amplitude of the oscillation). Output variables are 
a) , n and x , and B|_, the equivalent absorber admittance for the 
given geometry. 

Option 4 

This option is the same as Option 3 except that n and x are input 
and Bj is output. 

Option 5 

The last two options use B T as the iterated variable. They are most 
useful in generating stability maps in terms of n and x (Option 5) or 
real (Bj) and imag (Bj). (Option 6). Examples of such stability maps are 
presented in References 1, 2, 3 and 14. Frequency is used as parameter 
along these curves. 
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Option 5 is designed to compute Bj for a given complex frequency and 
a slot absorber. All the slot absorber parameters necessary for Option 3 
must be supplied here as well, in addition to the complex frequency. Output 
includes y^ n , Bj > n and x , and B L » the equivalent liner admittance. 

Option 6 

This option also uses Bj as the iterated variable. In this option, 
however, the absorber is characterized by an admittance, B|_ , and a length 
W a = X B - X A . For this option oj , S[_ , X^ and X-, must be supplied and 

output will give $j , y^ , and n and x . 

A summary of the principal input and output variables for the six 
options is given in Table 1 below. 


TABLE 1 


OPTION 

ONE 

TWO 

THREE 

FOUR 

OUTPUT 

1 

real (Bj) 

i ma g ( 3 j ) 

real (B L ) 

imag(3 L ) 

w , n & x 

2 

n 

T 

real (B L ) 

irrag(B L ) 

W , Bj 

3 

real (Bj) 

imag(Bj) 

BR 

AMF 

uj , " & x , 

4 

n 

X 

BR 

AMF 

w i 3 j » 3j_ 

5 

real (oj) 

imag(oj) 

BR 

AMF 

3 t , n & x , B 

6 

real (uj 

imag(od) 

real (B l ) 

imag(B L ) 

3j , n & x 


For convenience four main input variables are called ONE, TWO, THREE 
and FOUR both in the program and in the table. These variables represent 
different quantities in the different options. For example, in Option 1 
variable TWO represents the imaginary part of Bj , whereas in Option 2 it 
represents x , the time lag. 
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Choice of Fundamental Acoustic Mode 

As was mentioned in the "Theory" section, the success of the iteration 
process revolves around the assumption that the oscillation with active 
walls and mean flow is similar to one of the normal acoustic modes (no 
flow, hard walls) of the combustion chamber. The most useful variable for 
determining the suitability of an acoustic mode choice is the real part of 
the complex frequency. As a general rule, when the frequency nf oscilla- 
tion in the combustor is within 10% of a particular acoustic mode frequency, 
convergence will usually occur if that acoustic mode is used for 2^^ in 
the iterative process. Since the imaginary part of the frequency can be as 
large as the deviation of the frequency from its acoi'^tic value, nondimen- 

sional decay (or growth) rates as large as 0.20 (of the order of 1000 sec ^ 

★ # # 

for typical a* and R ) can occur for these conditions. 

For lower acoustic modes there is considerable separation in frequen- 
cies. At higher frequencies a given frequency may be close to two (or more) 
acoustic modes. In this latter case convergence problems can occur and it 
may be necessary to test all of the possible acoustic modes sequentially. 
Experience with the program must be the guide in these cases. 

For many (if not most) applications the acoustic mode choice is clear. 
For example, suppose that in a given combustor of diameter 2 ft, length 2 ft 
and average sound speed, i* , of 3000 ft/sec, an oscillation of frequency 
5700 sec”^ wore observed. The real part of the nondimensional frequency 
would be u R = (5700)/3000 = 1.90. This value is within 10% of 1.841, the 
acoustic frequency of the first transverse mode. (J^ O^r), l = 1, m = 1, 
n = 0, = A^) Hence, when investigating this oscillation using the 

iterative model, = ^ 110 (i.e., 


i-l,m«l,n»0). Indeed, the 
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choice of mode would be the same for frequencies between about 4970 sec”* 
and 6075 sec”*. On the other hand, if the observed frequency were 
6900 sec”* the oscillation would be closer in frequency to the combined 
first transyerse, first longitudinal mode frequency of 7260 sec”* (J-|(^i-| r ) 

cos T ’ 1 = 1 » m = *• n = '» n £mn = ( X ll + ij) ) and fl lll 

(£ = l,m = l,n = l) should be used. When the frequency of oscillation is 

"in between" two acoustic frequencies and is within 10% of neither, it is 

usually best to pick the higher mode. For the example given, if = 

6300 sec”* it would be within 10% of neither the pure first transverse 

frequency nor the combined first transverse, first longitudinal frequency. 

The best choice in this case would be rather than . 

The choice of acoustic mode is input to the program through the choice 
of the three integers, £ (radial), m (azimuthal), n (axial). In the 
program these are called I.HAT, MHAT, and NHAT, respectively. 
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Program MODULE 

In this section a general description of the program will be given 
first. Next, discussions of input and output formats will be given. 

Finally, a sample run will be presented and discussed, and a complete 
program listing will be given. 

General Description of Program 

The structure of the programs is as follows. (See Figure 3 for a flow 
diagram and Table 2 for a listing of program nomenclature.) After the non- 
default type variables have been declared, and the matrices and arrays have 
been dimensioned, the values for constants (such as PI fir]) are stored. 

Next, values for the two program variables K and IDCR are stated. The 
values for the iterated variable and the percent error in the modulus, and 
the absolute change in angle, will be printed out the first, last and every 
iteration. IDCR is an arbitrary number, after which a percent error in 
the modulus of over 50% will indicate that the problem is not converging. 

The values for the two constants K and IDCR (ID critical), may need to be 
changed, but it was felt they would not be changed often enough to warrant 
including them with the other input data. 


START 
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Flow Diagram 
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TABLE 2 

Computer Program Nomenclature 

Input Variables 
AMF - aperture mean flow 
ACAV 

BETAI - 
BETAL - 
BETAN - 
BR 


A c nondimensional sound speed in slot absorber cavity 


Bj, acoustic admittance of injector 

B l> acoustic admittance of liner 

3^, acoustic admittance of nozzle 


- backing distance 

EPSIL - e, amplitude of wave oscillation (only used in calculation 
of absorber resistance) 

ERROR - acceptable % error m magnitude and absolute change in 
radians between two successive iterations 

FOUR - OPTION = 1,2 or 6, imao(BE T U) 

OPTION = 3,4 or 5, AMF 

GAMMA - y, ratio of specific heats 

IDMAX - maximum number of iterations 

IN - n, interaction index 

K - iterated variable is output every K tn iteration 

LENGTH - L, nondimensional chamber length 

LHAT - $, radial acoustic mode, integer 

LTS - number of terms in radial direction 

MACH - M, Mach number 

MHAT - m, transverse acoustic mode, integer 
NHAT - n, longitudinal acoustic mode, integer 

NTS - number of terms in radial direction 
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ONE 

OPTION 

ROAP 

ROCAV 

TAU 

THREE 

TWO 

WCAV 

XA 

XB 

Program 

A1 

A2 

BES 

BETAIN 

CIOM 

DZP1L 

ETA 

ETA1 

ID 

LAMDA2 


- OPTION = 2 or 3, real (BETA I) 

OPTION = 2 or 4, IN 

OPTION = 5 or 6, real (OMEGA) 

- integer value between 1 and 6, sets which way the problem 
will be calculated 

- p^, density in slot absorber aperture 

- p c> density in slot absorber cavity 

- r , sensitive time lag 

- OPTION = 1, 2 or 6, real(BETAL) 

OPTION = 3, 4 or 5, BR 

- OPTION = 1 or 3, imag(BETAI) 

OPTION = 2 or 4, TAU 
OPTION = 5 or 6, imag(OMEGA) 

- nondimensional slot absorber cavity width 

- distance from injector to beginning of acoustic liner 

- distance from injector to end of acoustic liner. 

Variables 

- 4> evaluated at injector 

- (j> evaluated at nozzle 

. = /o j £(w ,rdr 

2X fcm 

- new BETAI 

- AM 

- || evaluated at liner midpoint 

• "L 

- n|* 

Zn 

- iteration counter 

- » here - ’ 0 
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f 


LEFF 

MU 

MUX 

NORM 

NORM! 

NPIL2 

OMEGAN 

PIL 

PIXL 

PS I 
PIL 
SINJ 

SLIN 

SNOZ 

VOL 

WA 

wwo 

XL 


y matrix 
old AMU 

in 

A 1 / 2 

in 

A 2 

(*£> 

new OMEGA 

TT 

L 

ttXL 

L 

y 

<f> evaluated at liner midpoint 

n dS iNj 

ii d \ IN 

II dS NOZ 

/// d V 

VI , width of liner aperture 

a 

OMEGA 

~ETKT 

distance from injector to midpoint of liner. 


Output Variables (not defined above) 

KI - Imaginary part of absorber impedance 

(reactance) 

RO - R real part of absorber impedance 
0 (resistance) 
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A point right after this (570 CONTINUE) is where the program returns 
to begin execution of a given set of data. The central processor (CP) 
time is stored at the beginning of execution of a set of data. This time 
is used to calculate the execution time for the set of data. Next, the 
data is read in and the counter ID is initialized. The constants for the 
set of data, such as (y-) , are then calculated. The first guess for the 
iterated variable and the p matrix (MU) are calculated, using the sepa- 
ration of variables solution. The input data and first guess are printed 
out. The setup is now complete, and each iteration returns to a point 
just below this (60 CONTINUE). 

For each iteration the iterated variable is first updated, then 
variables that are functions of the iterated variable are updated. The 
p matrix is stored in an extra matrix (MUX). The portion of the program 
from here to 500 CONTINUE, is designed to evaluate Equations (14) and (15), 
which give a new p matrix and iterated variable, respectively. The next 
section, down to 550 CONTINUE, does the following: Calculates the percent 

error in the modulus (ERR1), and the absolute error in the angle (ERR2). 
Then checks to see if the problem appears to be converging (ID greater than 
IDCR and ERR1 greater than 50%), or has gone the maximum number of itera- 
tions. If neither of the above has happened, the iterated variable, ERR1 

j. L. 

and ERR2 are printed out, if first or K tn iteration. Then checks to see if 
the problem has converged (ERR1 and ERR2 are both less than ERROR). If the 
problem has not converged, the program returns to 60 CONTINUE. If the 
problem is not converging, has gone the maximum number of iterations, or 
has converged, the final values for the iterated variable, ERR1 , ERR2, and 
the final p matrix, are printed out. If the problem does not converge, 
a message is printed out. 



The next segment of the program, down to 710 CONTINUE, calculates 
and prints the other information that is to be output. If a slot absorber 
is used, an equivalent BETAL is calculated. If IN (n) and TAU (t) are 
used, the equivalent BETAI is calculated; otherwise, the corresponding IN 
and TAU are calculated for the final OMEGA and BETAI. The CP time at the 
end of the program is stored. Running time is calculated and printed out. 
The calculations for this set of data are then complete and the program 
returns to 570 CONTINUE to begin calculations for the next set of data. 

If no more data is found, the program jumps to 6000 CONTINUE and stops 

> 

without an error message. 

The program is capable of handling up to 10 terms in the radial 
direction (LTS), 50 terms in the longitudinal direction (NTS), and a 
transverse mode as high as 4 (MHAT). This should be satisfactory for most 
cases. However, if the number of terms in the radial or transverse direc- 
tions must be increased, the Bessel values and Bessel roots for higher 
modes must be added to subroutines BESVL and BESRT, respectively. Also, 
all relevant dimension statements must be increased. To increase the 
number of terms in the longitudinal direction, only the dimensions of MU 
and MUX must be increased. If MHAT and NHAT are zero (0), LHAT cannot be 
one (1). This is a trivial case. 

The best compromise between good accuracy and fast running time (as 
discussed previously) occurs with 15-20 terms in the longitudinal direction 
(NTS), and 3 or 5 terms in the radial direction (LTS). LTS should always 
be odd, because the series has alternating signs in the radial direction. 

If the Mach number is greater than .40, more terms should be kept in the 
longitudinal direction. 
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Conmonly, for liners covering less than one-third of the chamber 
walls, evaluating the integral over the surface of the liner, by eval- 
uating at the midpoint and multiplying by the width, gives a good approx- 
imation to the integral with a big saving in running time. The program 
is set up to run this way. However, if it is desired to carry the inte- 
gration out, replace the SLIN card with the CALL LINER card and include 
SUBROUTINE LINER. (See program listing after 40 CONTINUE.) 

The final y matrix, which is printed out, should always be checked 
to be sure that the term corresponding to the acoustic mode assumed is the 
largest term in the matrix. If this is not the case, the wrong primary 
acoustic mode has been assumed, and the answer is not a characteristic 
of the primary mode assumed. 

Program Input 

The input necessary and the formats for typing the data cards are 
listed in Table 3. Before the first data card can be typed it must be 
decided from what is known about the engine (or desired from the calcu- 
lation) what value OPTION must take. It will be helpful to refer to the 
"Computational Methods" section in making this determination. Once the 
value of OPTION is fixed. Table 1 is used to determine which values the 
variables ONE, TWO, THREE and FOUR must take. The first data card can 
then be typed. 

The second data card contains the model information which must be 
supplied regardless of option. All inputs are real numbers. The third 
card contains program variable information, including convergence and 
matrix size limitation information. All variables on card three are of 
the integer type. The fourth data card needs to be included only when 
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table 3 


The first card: 

COLUMNS 

VARIABLE 

TYPE 

1-20 

ONE 

Real number 

21-40 

TWO 

Real number 

41 -60 

THREE 

Real number 

61-80 

FOUR 

Real number 

The second card: 

COLUMNS 

VARIABLE 

TYPE 

1-10 

real (BETAN) 

Real number 

11-20 

imag(BETAN) 

Real number 

21-30 

GAMMA 

Real number 

31-40 

MACH 

Real number 

41-50 

LENGTH 

Real number 

51-60 

XA 

Real number 

61-70 

XB 

Real number 

71-80 

ERROR 

Real number 

rhe third card: 

COLUMNS 

VARIABLE 

TYPE 

1-10 

LHAT 

Integer 

11-20 

MHAT 

Integer 

21-30 

NHAT 

Integer 

31-40 

LTS 

Integer 

41-50 

NTS 

Integer 

51 -60 

IDMAX 

Integer 

61 -70 

OPTION 

Integer 

he fourth card: 

COLUMNS 

VARIABLE 

TYPE 

1-10 

EPSIL 

Real 

11-20 

ROAP 

Real 

21-30 

ACAV 

Real 

31-40 

ROCAV 

Real 

41-50 

WCAV 

Real 

51-60 

LA 

Real 
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a slot type absorber is present in the combustor. All variables appearing 
on this card are real. 

When setting up the program the correspondence between text and 
program variables given in Table 2 will be useful. Also, it should be 
remembered that all variables in the program and text are nondimensional . 


Program Output 

The primary outputs of MODULE are the matrix y^ n and the iterated 
variable, either w or Bj • Values for these quantities are computed 
at every step and u (or Bj) is written out for every K^' iteration. K 
can be changed by the user by replacing a single card. In the program as 
presented here K = 5. The first and last iterations of y^ n are also 
printed out. The first iteration is a solution with no liner effect. 


consequently all terms in y^ n except y 


Jin 


are null entries. 


In addition, all model design variables are printed and labelled 
according to the program names of Table 2. The example to be discussed 
next demonstrates the typical form the output takes. 


Sample Run 

The combustor used for this sample run has a ratio of length to radius 
(LENGTH) of 2.0, a mean flow Mach number (MACH) of .3, and a ratio of 
specific heats (GAMMA) of 1.2. From the Mach number and ratio of specific 
heats, a nozzle response (BETAN) was calculated using the equation for a 
short nozzle given in the theory section of this report. The value for 
BETAN is .025 + O.Oi. There is an acoustic liner, of known admittance 
(BETAL), of .075 + O.Oi in place covering 10% of the cylindrical surface 
of the chamber, beginning one-tenth of the chamber length downstream of 



the injector face. Since space dimensions are nondimensionalized by 
dividing by the chamber radius, and the riondimensional chamber length 
is 2.0, this gives an Xa of .2 and an Xb of .4. The interaction index 
IN (n) and sensitive time lag TAU (t) are known to be equal to .30097 and 
2.219497, respectively, for the injector •'sponse. With the information 
known about the injector and liner responses and looking at Table 1, it is 
determiner* that OPTION must equal 2, and ONE is IN, TWO is TAU, THREE is 

real (BETAL), and FOUR is imag (BETAL). Table 1 shows that (w; OMEGA, the 

nondimensional frequency and decay rate, and the effective FFTAI for the 
injector will be calculated. The first transverse acoustic mode is chosen 
as the primary mode. This corresponds to LHAT = 1, MHAT = 1, and NHAT - 0. 
A good compromise between running time and accuracy was desired, so by 
referring to the discussion in this report, the number of terms chosen in 
the radial direction (LTS) is equal to 3, and the number of terms chosen 
in the longitudinal direction (NTS) is equal to 16. A high precision is 
desired, so ERROR is chosen as .01. The maximum number of iterations 
allowed for this set of data (IDMAX) will be 50. Since an acoustic liner 
is used, no fourth data card is needed. The input values are then typed 
up on three data cards, as described in Table 3. Tm cards as punched and 
submitted appear in Figure 4. 

The output from the sample run is shown in Figure 5. The pregram 
first prints out the value of every variable that is input on the data 
cards. This is to allow for double checking, to be sure all the input 
data is correct, and also so there is a complete description of the rocket 
engine that was simulated. Next, the fundamental frequency for the primary 

mode assumed, and the first guess of the iterated variable, and tue u 
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matrix are printed. (These are the separation of variables values.) If 
OPTION is 5 or 6, the first guess of BETAI is printed under the injector 
response description. Next, the value of the iterated variable is printed 

i. u 

for the first, every K , and last iteration. In the sample run K =5, 
so the first, fifth and sixth (last) iteration values are printed, along 
with the errors in the modulus and phase angle of the iterated variable 
(w in this case). The last matrix (iteration 6) is then printed 

out. It can be seen that the term corresponding to l = i = 1 and n = n = 0 
is, indeed, the largest. In fact, in this case it is an order of magnitude 
larger than any other matrix element. The other output information is then 
printed. In this case, t.;e BETAI calculated from the input IN and TAU and 
the final OMEGA. The last thing printed out for each set of data is the 
begit.ning time TBG, ending time TEND, and execution time TEX, for this set 
of data. 

Program Listing 

A complete program listing is presented at the end of this report. 
Comment cards are used liberally and much of the program is self-explanatory. 
The computer program MODULE conforms to Fortran IV ANSI standards. 


t 
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Figure 4. Data Cards for Sample Run 
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THF LENGTH 

Trt RADIUS = 

P.000000 

hF an FLUr hach NO. * .300000 

LHAT s 1 

mhat = 1 

NHA T = 0 

(iPTIOH = t 

■SA*M A = 1 . 

pftuoon 




THE ABSOPPF 0 HAS THE FOLLOWING GEOMtltf 

th^ lil*-w starts »t .pooooo and fni>s at .4onn<m ap*>thpf width s .po^ono 
F-i-TAL = .07^000 U.OHOOCO 
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0. 

II • 

(I. 


0 

M 

0 

0 

0 

<1 

n 

II 

N 

n 
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Figure 5. 


Output for Sample Run 
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ANf » l _* 


i 

* 

6 


1.799*4 
1.78*49 
1 .78837 


.18237SE-01 

• 1777A7F-m 

• 1 767?v*»i-0 1 


1.0077** 
.034 1 *<- 
. 0f»8“>33 


.01 013? 
. or . o ? 3 S 
. 0000*1 


THF FINAL M U MATRIX JTS AS FOLLOWS 


N 


L = 1 


0 

1 

? 

3 

4 

5 

6 
7 
h 
Q 


10 

11 

IS 

IS 


2.04*8? 

.2*3107 

.?*39**E-0l 

. 70677SE-U? 
.44931*E-u2 
• 3SSS32E-0? 
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-. Mul * 9 c - O ? 

- • *Sllft J 1 K - 0 /» 
-.21?9*3^-n; 
-.331 1'^7^ - 0 ^ 
- • ] 94 * 7 rt - - n * 

- . ?9t),j7**- -•/? 

- • P f > — ** ^ 

-.?* t '^02F-o? 
-. J 7*8OHF-02 

194**76-0? 
1 2 l 3 84 K — 02 


. 1 77y43F-03 

.7MrbU-03 
. M * JlM ' i-OJ 
. 3 i » Sr >^« F -0 3 
. 4 *»»* o 7 K -04 
-. 1 1 /743F-03 

-.179 482K-0 j 
1 714*hh-<»3 

- • ?06n0 ,£-04 

- • 1 * h ^ 40 F - l *4 
.223S23F-04 
.4**V v9t-04 
. S 077 '» OL -04 
.474f> J h£-()4 


L=2 


.S334ift»--ii? 

.874*4 V-o? 

• 4b*7 3?r — 0? 
.94030^*’ -03 

-.13*034*-'.? 
-.23?* 32F-0? 
-.23901 7 P - » V 
-.1434*3^-.'? 

- » l ? ft 4 W 4 *> - U ? 
-.S720K Hr -03 
. 1138 V , Hr -»:4 
.4l*l?4- -l»3 
.^ 24 H ^ hF -(»3 

• * S 4 . 0 ft 0 r - n 3 

• *47fiftOi» -03 

• Jd * 904 J - - u 3 


-.1*41 -I Sf.-ni 

-. 31*7141- -03 

-.?*3?v fr -0 j 

-,??l4^ftf -03 
- . 44 / 7 >41- -I|4 
.29 04 
. 1 ?•* 7M*- - 03 

. 1 -03 

. 1 1 , >*n 0 »)F — ( 1 3 

.809?73K -04 

. 1 S*ft 1 -(14 

2?**X HI -04 
-.*♦ H 09 l *< F -04 
-.9?2^9?f-r.4 


1=3 


?3Mi>*-02 
40M9^7 f - 02 
24ft7j7h-0? 
‘'> 7 h 84 i . i - - u 3 
47ovjft^i -03 
l ' H * n ? fiF - ii ? 
20M O 7«h-0? 
1*1) 74,. - 0 ? 
1 ? 4 S 8 *»- - II ? 
ftMblSftr -(13 
ft 4 Ort??F — O * 
44 - 03 

^9?4 * J c - f* 3 
72ft77?F-ij 1 
HlShoSi^ -03 
4(U73SF-03 


The FIN4L RFTaI IS .o<xM* 1 .0iSR17h 

TfiCis H.A10 TFM) = IS.S 45 TFxs !.)<;? 


95**MI. PAQC 13 

orp °msmm 


Figure 5 


( Continued) 


- 40 - 


<i 4 


C 

C* 

c« 

c* 

c* 

c» 

C* 

c* 

c* 

C* 

c* 

C* 

c* 

c* 

c* 

c° 

c* 

c* 

c* 

C* 

c« 

C* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c» 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c* 

c« 

c* 

c* 

C* 

c* 

c* 


PROGRAM MOOULF ( INPUT * OUTPUT * T APf 5= INPUT * T APT 6-OUTPUT) 
>>SM>'0<M>0O646att0<M»Oa00ee0e6aOaa00#aa404}ae6&0Oe0eO00ee0e00iK0tf<}{>&4HMH>44ee 

COMPU.'EH PROGRAM MOOULF 

WRITTEN AT COLORADO STATE UNIVERSITY 

DEPARTMENT OF MECHANICAL ENGINEERING 

f-ORT COLL INS » COLORAOO ByS23 

SPONSORED 8Y NASA LEWIS RESEARCH CENTER 

GRANT NGR 06 - 002 - 09S 

DIRECTED BY RICHARD j. PRIEM 

™IS PBOORAM IS DOCUMENTED IN THE MASTERS THESIS OF KURTIS W. ECKERT 
l?I£A.£2k 0HAL)0 STArE UNIVERSITY* AND CONFORMS TO ALL FORTRAN IV ANSI 
b I ANOARD5 • 

IT IS WRITTEN TO GIVE A LINEAR ANALYSIS OF HIGH FREQUENCY 
COMBUSTION STABILITY IN LIQUID PROPELLANT ROCKET ENGINES. THE 
P^ Y ^ICAL MODEL USED IS A RIGHT CIRCULAR CYLINDER. THE COMBUSTION 
IS MODELED AS EITHER AN ARBITRARY ACOUSTIC ADMITTANCE (BETAl). OR 
HT JHE CROCCO SENSITIVE TIME LAG THEORY (IN AND TAU). THE N02ZLF 
}§ W^OELEO AS AN ARBITRARY ACOUSTIC ADMITTANCE (BETAN). THF LINER 
IS MODELED AS *N ARBITRARY ACOUSTIC ADMTT4NCE (HETAL> OvfR 

5>OHf PORTION OF THF CYLINDRICAL WALL OF THE CHAMBER. OR AS A SLOT 
ftBSOHBEri • 

this rrogram has rpen writtfn to solve FOR THF nondimfnsional 
complex frequency if all the boundary responses arf given. OR solve 
for the INJECTOR WESPONSE IF The COMPLEX FREQUENCY AMD L I NFR AND 

§I VtN » THIS LEADS TO A TOTAL OF 6 OPTIONS 
FOR RUNNING THE COMPUTER PROGRAM. 

Tim , POLLOWING IS A TABLE TO USE IN DETERMINING WHICH VALUE OF OP- 
TION TO USE. THE TABLE SHOWS WHAT INFORMATION MUST BE GIVEN AND 

CALCULATED FOR EACH OPTION. R(> MtANS THE REAL PART OF 
THE COMPLEX VALUE INSIDE THE PARENTHESES. AND Id MFANS THE IMAGIN- 
ARY PART OF THE COMPLEX VALUE INSIDE THE PARFNTHESES. 


TABLE OF OPTIONS 
OPTION 


INPUT VARIABLES 
ONE TWO THREF 


1 

3 

4 

5 
b 


R ( BETA I) 
IN 

R * BET A 1 ) 
R (OMEGA ) 

r (omega; 


I (BFTAI > 
tau 

I (BFTAI ) 
Tau 

I (OMEGA) 
I (OMEGA) 


R(BETAL) 
R (BETAL) 
HR 
BP 
BP 

R(BETAL) 


FOUR 

1 (BETAL) 
I (BETAL) 
AMF 
AMF 
AMF 

I (BETAL) 


OUTPUT VARIABLES 


OMEGA. IN & TAU 
OMEGA. BETAl 
OMEGA. RFTAL. IN E TAU 
OMEGA. BFTAI. BETAL 
BFTAI. 0FT4L « IN N TAU 
BETA!. IN L TAU 


AK. «r»F2 L k9rI N 9 Al A LISTING THE INPUT VARIABLES THAT ARE PUT 

ON EACH DATA CARD. AFTER THE VARIABLE NAME THE TYPE OF VARIABLE IS 
r22 W r.*n T AF N r A r- BW JF F DESCRIPTION OF THE VARIABLF IS GIVEN, THEN AT 
THE END OF THE LINE ARE THE COLUMNS OE THE DATA CARD WhICH THE VALUE 
FOR THIS VARIABLE MUST BE TYPED IN. ALWAYS BE SURE TO RIGHT JUSTIFY 
INTEGER VALUES. THE FOURTH DATA CARO IS USED ONLY WHEN OPTION EQUALS 
3. 4 OR 5. 


4 

p 

LIST OF 

INPUT 

4 

F IRST CARO 

4 

ONE 

REAL 

o 

TWO 

REAL 

4 

three 

REAL 

4 

P 

FOUR 

REAL 

4 

SECOND 

CARD 

4 

4 

BETAN 

COMPI 

4 

4 

GAMMA 

REAL 


SEE TABLE ABOVE 
SFE TABLE AROVE 
SEE TABLE ABOVE 
SEF TABLE ABOVE 


ACOUSTIC ADMITTANCE OE NOZZLE 
MODULUS LESS THAN .5 
IF NO KNOWN VALUE USE SHORT NOZZLE 
RATIO OF SPECIFIC HEATS 


1-20 

21-40 

41-60 

61-80 


1-10 6 , 11-20 
21-30 


t 
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c* 

MACH 

REAL 

MF AN FLOW MACH NUMBER 0 < MACH < .S 

31-40 

c* 

LENGTH 

REAL 

LENGTH OF CHAMpER/RADIUS OF CHAMHF R 

A I -50 

c» 

XA 

REAL 

DTSTAiJCE from injfctor face TO start 

51-60 

C* 



OF LlNfc'R/RAOIUS OF CHAMBER 

c* 



0 < XA < XR 


c* 

XH 

REAL 

DISTANCE FROM INJECTOR FACE TO END 

61-70 

c* 



OF |. INER/PADIUS OF CHAMHFR 

c* 



XA < XB < LENGTH 


c* 

c« 

c* 

c* 

c* 

c* 

FPPOR 

real 

MAXIMUM ALLOWABLE % ERROR IN MODULUS 
OR ABSOLUTE DIFFEPANCE IN RADIANS OF 
ANGLE TO DETERMINE CONVERGENCE OF 
ITERATED VARIAHLF 10E-5 < ERROR < 1 

71-80 


THIRD CARD 



c* 

LHAT 

integer 

ASSUMED MODE IN RADIAL DIRECTION 

1-10 

c* 



TYPICALLY 1 OR 2 

c* 

MHAT 

INTEGER 

ASSUMED MODE IN TRANSVERSE DIRECTION 

11-30 

c* 



TYPICALLY 0. 1 OR 2 

c* 

nhat 

INTEGER 

ASSUMED MODE IN LONGITUDINAL DIRECTION 

21-30 

c* 



TYPICALLY 0. 1 OR ? 


c* 

LTS 

INTEGER 

NUMRFR OF TERMS KEPT IN RAOIAL DIRECTION 

31 -AO 

c* 



ODD < 10 TYPICALLY 3 OR 5 


c* 

NTS 

integer 

NUMBER of terms kept in long, direction 

<*1-50 

c» 



NTS < SO TYPICALLY IS < NTS < 20 


c* 

IDMAX 

INTEGER 

MAXIMUM number of iterations allowfo 

51-60 

c* 



TYPICALLY 100 

C* 

c* 

c* 

c« 

OPTION 

INTEGFR 

SFE table, above 

1. 2* 3, A. S, OR 6 

61-70 

fourth 

CARD 



c* 

EPSIL 

REAL 

WAVE AMRL T TIJOE 

1-10 

c* 

ROAP 

RtAL 

APERTURE DENSITY ratio 

1 1-20 

c* 

AC A V 

PEAL 

CAVITY SOUND SHE Eli 

21-30 

C* 

ROCAV 

real 

CAVITY OtNSITY RATIO 

31-40 

c* 

WCAV 

REAL 

CAVITY WIDTH 

A 1 -50 

c« 

LA 

REAL 

APFRTURE LFNGTH 

51-60 


(.» 

c 

REAL INI . I NO, IN,NPIL?.LE.NGTH,LEEE,LA,L AMDA.L AMDA2.N0rM.maCH,mACH2 
1 • NORM 1 

INTEGER OPTION. CHECK 
C 

COMPLEX HETAN.BE TA| . BE T A I , MU . MIJX .OMEGA ,UMfc GAN. C I OM.CFRROR , C I 
1 .CZfcHO.F 1 .HLINEK.Pit .DZPll .PRES.Al « A2.SUM1 .SUM?. 

\ ?AV. BASF*. RETAIN, WWO 

% COMPLEX CTFRM.H1.B?.H1SU.H?^0.A,FXP1 ,EXP?.TFRM1 .TFRMP.VO! .S1NJ, 

1SNOZ, SINJB.SLIN.PS! 

C 

DIMENSION MU ( SO. 10) .MUX ISO. 10) .A1 I 111) . A?( 10) 

INITIALIZE CONSTANTS 

PI =3. 1415926535 1 ) 

CI=CMPL* 10. . 1 . ) 

C7ERO = CMPL X ( 0 . . 0 . ) 

I0Ch= ?0 
K=5 

PF AO IN DATA 

B70 CONTINUE 

CALL SECOND I TOON ) 

READ (S.100) ONE .Two, T hree .FOUR. HE T AN. GAMMA .MACH.LE NOTH. X A . AH • 
lFRROR.LHAT.MHAT.NriAT.LTS.NTS, I UMAX .OPT ION 
IF(EOF(b)) BO00.S80.6n0O 
580 CONTINUE 

IF (OPTION. GF. 3. AND. OPTION. LE. 5) READ (5.101) EPSIL .ROAD, 

1 ACAV.RUCAV.WCAV.LA 

INITIALISE VARIABLES 






A 
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ANG1 = 0.0 
ID = 0 

IF (OPTION. E0.1. OR. OPTION. EQ. 3) 8ETAIN= CMPt X (ONE , TWO ) 
I Ns ONE 
TAU= TWO 

RETAL= CMPLX(THREE.FOUR) 

BP= THREE 
AMF= FOUR 

CALCULATE CONSTANTS 

P IL= PI/LENGTH 

NPIL2= (FLOAT (NHAT) *PIL> 

WA = XB -XA 
XL = XA ♦ WA/2.0 
PIXL=PI*XL/LENGTH 
LAMDA2=BESRT (MHAT.LHAT) **? 

IF (OPTION. LE. 2. OR. OPTION. EO. 6) GO TO 30 
LEFF = LA ♦ O.37B*O.05*WA*( 1. - 0.7 »SQRT< WA/WCAV) ) 
30 CONTINUE 

ETA1= SORT (LAMDA2 ♦ NPIL2) 

CALCULATE FIRST GUFSS FOR ITERATED VARIABLES 

OMEGANs .98«CMPLX (ETA1 .0.0) 

IF (OPTION. GE. 5) 0MEGAN= CMPI.X (ONE. TWO) 

W*0=0MEGAN/ETA1 

IE (OPTION. EQ. 2. OR. OPTION. EQ. A) BETAIN=MACH# ( 1 ./GAMMA 
1 E XP ( -C I ^OMEGAN^T AU) ) ) 

CHECK VALUES OF INPUT VARIABLES 


TERM= WA^REAL (PETAL ) 

CHECK= 0 

IF (MACH.LE.O.O.OR.MACh.GE.I.O) CHFCK = 1 
IF (MACH.GT.0.S) WRITE (6.Q00) MACh 
IF (GAMMA. LT. 1 .0. OR. GAMMA. GT. 1.67) CHECKS l 
IF (GAMMA. LT. 1 . 1 ) WRITE (6.920) Gamma 
IF (LENGTH. LE. 0.0) ChECK= 1 
IF (LENGTH. GE. 3.0) WRITE (6.902) LENGTm 
IF (LENGTH. LE. 1.0) WRITE (6.92?) LENGTH 
IF (REAL (BETAN) .LT. 0.0) CHECKS 1 
IF (REAL (BETAN) ,GT. 0.3) WRITE (6.906) HtTAN 
IF (OPTION. GT. 2. AND. OPTION. LT. 6) GO TO 44 
IF (TERM. LT. 0.0) CHECKS I 
IF (TERM. GT. 0.3) WRITE (6.904) TERM 
44 CONTINUE 

IF (OPTION. NE. 2. ANO. OPTION. NE. 4) GO TO 46 
IF (IN.lT.O.O.OR.TAU.LT.O.O) CHECK= 1 
IF (IN.GT.J.O) WRITE (6.91?) IN 

IF (TAU.GT.4.0) WRITF (6.9J4) TAU 

46 CONTINUE 

IF (OPTION. NE. 1 .ANO. OPT ION. NF. 3) Go TO 47 
IF (CABS (RET A IN) .6T.2.0) WRITE (6.916) HFTaJN 

47 CONTINUE 

IF (XA. LT. 0.0. OR. XA.GT.XH. OR. XB.GT. LENGTH) CHFCK- \ 

IF (WA.GT.O.F) WRITE (6.93(1) WA 

IF (LHAT .LE. 7. ANO.MHAT . LE < 4 . AND.NHA T . LF . 1 0 ) GO TO 4 H 
CHECKS 1 

WRITE (6*980) LHAT .mmaT.NHAT 

48 CONTINUE 

IF (NTS.LT.NHAT ♦ 10) WRITE *6,940) NTS.NHAT 
IF (LTS.LT.LHAT ♦ ?) WRITE (6.942) LTS.LHAT 

IF (REAL (WWO) .GT. 1 . 1 ) WRITF (6.9S0) 

IF (OPTION. I. T. 3. OR. OPTION. FO. 6) GO TO 49 
IF (EPSIL.LT.O.O.OP.HR.LF.O.O.OR.AMP.LT.O.O) CHFCK r ) 

IF (PR.GT.O.?) WRITE (6.908) HR 

IF (AMF.GT.0.2) WRITE (6.910) AMf 

49 CONTINUE 

IF (CHECK. EO. I ) WRITE (6,960) MACH, GAMMA .IFNGTH, HF T AL . RF T AN . XA , XR , 
1 I N . T A U 

IF (CHECK. £0.1) GO TO 710 


4 
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50 


INITIALIZE FIRST GUESS OF mu MATIX 

00 50 N1=1«NTS 

no so l=i*lts 
MU<N1.L>= CZEPO 
CONTINUE- 


70 


75 

71 


THF TILDA SOLUTION 

BV = 8ESVL (MHAT.l HAT) 

EP= 1, 

IF (NHAT.EQ.0) EP= 2 . 

RES1= RV*BV*(LAMDA2 z FtOAT (MHAT) **?) /2./LAMOA? 

NORM 1 = SORT (BES1 *EP*LEN6TH/?V) — __ 

MACH2= MACH*MACH 
DO 75 N= 1,100 

C TERM= MACH2*0MEGAN*0ME6AN* (MACH2-1 . ) * (L AMf^-OMFGAN^OMFMN ) 
CTERM= CSQRT(CTERM) 

B 1 = (MACH*OMFGAN»CTFPM) / ( 1 .-MACH2) 

82= (MACH*OMFGAf\'-CTFPM) / ( 1 ,-MACH2) 

R1S0=B1*B1 

B2SQ=B2*82 

A=- (CEXP(CI»Lt.NGTH* (R1-B2) ) * (HI ♦RF T AN*GAMMa* (OMFGAN»MACH»R 1 ) ) / 

I (P2*BETAN*GAMMA* (OMFGAN*r;ArH*R?) ) > 

EXPl=CEXP(CI«LENGTH*fll ) 

EXP2=CEXP(CI*LfcNGTH*B2> 

TERM1= B1*(EXP1*(-1.) **NHAT-1 . ) / (HJSQ-NPIL?) 

TEPM2= A*B2* (EXP?* (-1 . ) **NHAT-1 . ) / (RPSQ-NPIl 2) 

VOL= MACH* (MACH*(B1S0*TEPM1 *R?SO*TFPm?) ♦ 2.«0MFGAN* (RlttTFWMl ♦ 
1B2*TERM2) ) 

S I NJ=GAMMA* (OMEGAN* ( 1 .♦A) *MACH* (Hi »A»( ?) ) 

SINJ8=BETAIN*SINJ 

SNOZs BETAN*GAMMA » (-1 . ) **NHAT* (OMF GAN* ( fc'XPl ♦ A *t XP2 ) *MACH* ( R 1 * 
1EXP1 ♦ A*B2*EXP?) ) 

IF (OPTION, GE. 5) GO TO 70 

OMEGAs -(VOL ♦ SINJH ♦ SN07)/(TFWM1 ♦ TERM?) ♦ FTA1**2 
OMEGAs CSQRT (OMEGA) 

IF (OPTION. EQ.2. OR. OPTION. FO, 4) RETAIN= MAC W * ( 1 ./GAMMA-IN* ( 1 .-C 
1EXP (-CI*0ME6AN*TAU) ) ) 

CERRORs OMEGA - OMFGAN 
OMEGANs OMEGA 

IF ( ABS(REAL(CF.«ROP) ) .GT. 0.0001) 60 TO 75 
IF ( ARS(AIMAG(CERROR) ) .GT. 0.0001 ) GO TO 75 
GO TO 71 
CONTINUE 

RETAINS ( (ETA 1 **2- OMFGAN**?) *< 7 ERM1 *TEPM2> - VOL - SNOZ)/SINJ 

60 TO 71 

CONTINUE 

CONTINUE 

PSI= ( -2. *C I /LENGTH) * ( TEWM1 ♦ TFWM2) 

IF (NHAT.EQ.0) PSI= PSI/2. 

MU (NHAT* 1 ,LHAT ) = CMPLX ( 1 . /NORM1 ,0.0) 

DO 55 Nisi. NTS 

Ns N 1 — 1 
RNs FLOAT (N) 

IF (N.EQ.NHAT) GO TO 55 
TERMs FLOAT ( (-1 ) **N) 

NPIL2= (RN*PIL)**2 
FP= 1. 

IF (N.FQ.O) EP= 2. 

NORM= BES1*FP*LE NGTH/2. 

ETAs LAMDA2 ♦ NPIL? 

TFRM1= (EXP1*TERM - 1.)/(B1SQ - NPIL2) 

TERM2= (EXP2*TERM - l.)/(B25Q - NPIL2) 

VOL= MACH2* (fllStl*Bl*TERMl ♦ A*B2SQ*B2*TER -2) ♦ 2.*MACH*0MF6AN*(R 

I I SQ*TERM1 ♦ A»B2SO*TERM2) 

SN0Z= BETAN*GAMMA*TERM* (OMFGAN* (EXP1 ♦A*EXP?) ♦ MACH* (R 1 *E XP 1 ♦ 

1 A*B2*EXP2) ) 

SINJB= BETAIN*GAMMA*(0MEGAN*(1. *A) ♦ MACH*(B1 ♦ A*B2)) 
l/NORM/PSI T> * 8ES1 * CI * <V0l * SN0Z ♦SINJB)/(0MEGAN**2 -ETA)/N0HM1 
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55 CONTINUE 

PRINT OUT PPOBLEM DE SCP I PT T ON AND INITIAL VALUES 

WRITE (6*1) 

WRITE (6.200) 

WRITE (6.210) LENGTH. MACH. LHAT.MHATtNHAT.OPTION, GAMMA 
WRITE (6.220) 

WRITE (6.230) AA.XR.WA 

IF (OPTION. GE. 3. ANO . OPT I ON.LE .5) GO TO 20 
WRITE (6.235) BETAL 
GO TO 22 
20 CONTINUE 

WRITE (6.240) BR.AMF.LA.LEFF.EPSIL.WCAV.ROCAV.POAP.ACAV 
22 CONTINUE 

WRITE (6.250) 

IF (OPTION. GE. 5) GO TO 10 
IF (OPTION. EQ. 2. OR. OPTION. FO. 4) GO TO 12 
WRITE (6.260) RETAIN 
GO TO 14 
12 CONTINUE 

WRITE (6.270) IN, TAU. RETAIN 
GO TO 1 4 
10 CONTINUE 

WPITE (6,280) OMEGAN.WWO 
WRITE (6.285) RETAIN 
14 CONTINUE 

WRITE (6,290) 

WRITE (6.310) BETAN 

WRITE (6.320) LTS.NTS.IDMAX.ERROR 

WRITE (6.300) ETA1 

WRITE (6.305) OMEGAN 

WRITE (6,105) 

00 35 Nsl.NTS 
NM 1 s N - 1 

WRITE (6.104) NM1 , (MlMN.L) .L=l .LTS) 

35 CONTINUE 
WRITE (6,2) 

IF (OPTION. LF. A) WRITE (6.350) 

IF (OPTION. GE. 5) WPITF (6,360) 

BEGIN ITERATION 

60 CONTINUE 

UPDATE ITERATED VARIABLES ! 

IF (OPTION. EO. 2. OR. OPTION. EO. 4) HETAIN=MACh* ( I ./GAMMA - IN«>(1. -C 
1EXP(-CI*0MEGAN*TAU) ) ) 

BETA I = BETAIN 
OMF GA=0MEGAN 
CIOMsCI*OMEGA 
WR= REAL (OMEGA) 

STORE NEW MU MATRIX IN EXTRA MATRIX 

00 8u0 Nisi, NTS 
DO 800 L=1.LTS 
MUX(N1 ,L)sMU(N1,L) 

800 CONTINUE 

CALCULATE PHI AT INJECTOR. N07ZLE .MIDPOINT OF LINER 

R 1L= CZERO 
DZP1L* CZERO 
DO 130 L S 1 »L TS 
SUM 1 = CZERO 
SUM2= CZERO 
BV= BESVL (MHAT.L) 


C 

C 
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00 120 N1=1,NTS 
N= N1 - 1 
RN= FLOAT (NJ 
AV=MU(N1 »l> 

SUM1 = SUM1 ♦ A V 

SUM2 = SUM2 ♦ AV*FL0AT ( (-1 ) **N) 

P1L=P1L*BV*AV*C0S(RN*PIXL> 

DZP1L= DZP1L ♦ HV*AV*RN*PIl*SIN(RN*PIXL> 

120 CONTINUE 

A1 (L)=SUM1 
A2(L)=SUM2 
130 CONTINUE 

CALCULATE BLINER FOR MIDPOINT OF LINER 

PRES= GAMMA*(CIOM*PlL - MACH*07P1L) 

BLINER = WA*PRES*BFTAL 

CALCULATE BLINER FROM LINER GEOMETRY 

IF (OPTION. LE. 2. OR. OPTION. FQ. 6) GO TO 3000 

AKl = leff*wr*rgap-wa*acav*pocav/wr/wcav/hr *acav 
PRE1 = CABS (PRES) 

RO=l. 

ROl = RO 

DO 133 1=1.100 

F = SORT ( 1.0 ♦ (AK1/RO)**? ) 

BASE= EPSIL°PREI/F /GAMMA 

RO = SORT ( 0 . 8* ( 1.5»AMF*RO ♦ BASF ) ) 

R02 = ABS < (ROl-RO) /ROl ) 

IF(R02.LT.1.0fc-04) GO TO 134 
RO 1 = RO 

133 CONTINUE 

134 CONTINUE 

BLINER = WA*PRES/GAMMA/ <RO«CI»AKl ) 

3000 CONTINUE 

IF (OPTION. LF. 4) GO TO 3001 

CALCULATE NEW HETAT 

VOL = CZERO 
00 45 N1=2.NTS 
NM1= Nl-1 

IF (NM1.EQ.NHAT) GO TO 45 

Kl* NM1 ♦ NHAT 

K2= NMl - NHAT 

L1=(-1)**K1-1 

L2=(-1)«*K2-1 

LSUM=L 1 *L2 

IF (LSUM.EQ.O) GO TO 45 
Cl= FL0AT(L1>/FL0AT(K1> 

C2= FLOAT (L2)/FL0AT (K2> 

SC= (Cl ♦ C2) *f LOAT (NMl ) 

VOL= VOL ♦ MUX (N1 .1 NAT) "SC 
45 CONTINUE 

VOL= V0L < *MACH«CI0M*»FS1 

V0L= VOL - MUX(NHAT*l.LHAT)'>Ht51'»(NAC.N<»t-LOAT(NHAT)*PlL> ,,<> 2*LFNGTH 

1 / 2 . 

VOL= V0L/N0PM1 

SN02= BETAN*GAMMA4CI0M*BES1*A2(LNAT) ***LOAT ( ( - 1 ) *4NmAT ) / NORM 1 
SL IN= BL INER*BESVL (MmaT .LHAT ) *C0S (FLO* T (NHAT ) *PIXL ) /NORM 1 
SINJ= GAMMA*CIOM*0FS1*A1 (LHAT) /NORM1 
BFTAIN= (OMFGA**2 -FTA1**2 -VOL -SNOZ -SLIM/SINJ 

3001 CONTINUE 

START DO LOOP FOP L SUMMATION 
00 S00 L=1.LTS 
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« 

» 

« 


8 

C 


LAMDA* 8ESRT (MHAT »L ) 

BV2* BESVL(MHAT*L>**2 

START 00 LOOP FOR N SUMMATION 

DO 500 NX=1 *NTS 

INITIALIZE AND CALCULATE CONSTANTS FOR THIS SUMMATION 

N* NX - 1 
VOL* CZERO 

ETA* LAMDA ♦ (FLOAT (N) *PIL ) **2 

RES* BV2* (LAMDA - FLOAT (MHAT) **2> *.5/LAMDA 

FP* 0,5 

IF (N.EQ.0) EP= 1. 

NORM= BES*EP*LEN6TH 
N0RM= SORT (NORM) 

BASEW*OMEGA*OMEGA-ETA 

CALCULATE PROPAGATION TERMS 

DO 40 Nls 2, NTS 
NM l = Ni - 1 

IF (NM1.EG.N) GO TO 40 

K 1 = NM1 * N 

K2* NM 1 - N 

Ll*(-l)**KI-l 

L2-(-l)»*K2-l 

LSUMsL1*L2 

IF(LSUM.EQ.O) GO TO 40 
Cl* FLOAT (Li ) /FLOAT ( K 1 ) 

C2= FLOAT (L2) /FLOAT (K2> 

SC= (Cl * C2) *FL0AT (NMD 
VOL* vOL ♦ MUX ( N 1 * L ) *SC 
40 CONTINUE 

VOL* VOL*MACH»CIOM*BES 

VOL* VOL - MUX(NX»L)*BES*(MACH*FLOAT(N)*PlL) -»*2«LENGTH/2, 

VOL* VOL/NORM 

CALCULATE NOZZLE INTEGRAL 

SNOZ* BETAN*GAMMA*C10M*BES*A2 (L) *FLOAT ( (-1 ) **N) /NORM 
CALCULATE LINER INTEGRAL 

TO CARRY OUT THE INTEGRATION OVER THE LINER USE THE FOLLOWING TWO 
CARDS. AND PUT A »C" IN THE FIRST COLUMN OF THE CARD FOR APPROX- 
IMATING THE INTEGRAL * HENCE MAKING IT A COMMENT CARO. 

CALL L INER ( XA. XB.PIL.L.N. MHAT. MUX *CIOM.BLlNER.LTS»NTS*CZFRO) 

SL INs GAMMA*BESVL (MHAT .L > *BET AL*BL INER/NORM 

TO APPROXIMATE THE LINER INTEGRAL BY EVALUATING AT THE MIDPOINT 
AND MULTIPLYING BY THE WIDTH* USE THE FOLLOWING CARD. AND PUT A 
"C" IN THE FIRST COLUMNS OF THE TWO PRECEDING CARDS FOR CARRYING 
OUT THE INTEGRATION* HENCE MAKING THEM COMMENT CAPOS. "SUBROUTINE 
LINER" NEEO NOT BE COMPILED IF THE INTEGRAL IS TO BE APPROXIMATED. 

SL IN* BL INER*BESVL ( MHAT *L ) *COS (FLOAT (N) *P I XL > /NORM 

CALCULATE INJECTOR INTEGRAL 

SINJ* GAMMA*CI0M*BES*A1 (L)/NORM 
SINJB* SINJ*B£TAIN 

CALCULATE NEW MU TERM 


5F 


VOL ♦ SINJB ♦ SLIN ♦ SNOZ 
(N.EO.NHAT.AND.L.EO.LHAT) GO TO 430 


MU (NX *L > * Fl/NORM/BASEW 
GO TO 500 
430 CONTINUE 


k 
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C 

C 

C 


FOR PRIMARY MODE CALCULATE MU TERM AND NEW OMfGA 


C 

C 

C 


C 

C 

C 


MU (NX »L ) * CMPL X ( 1 . /NORM* 0 . 0 ) 

IF (OPTION. LE. 4) OMEGAN= CSQRT(F). ♦ ETA) 

500 CONTINUE 
ID= ID ♦ 1 

IKs IFlX(FLOAT(ID>/FLOAT(K)> 

RKs FLOAT (ID)/ FLOAT (K) 

RK= RK - FLOAT (IK’ 

IF (OPTION. LE. 4) GO TO 4000 

CHECK FOR CONVERGENCE ON BETAI 

CERRORs BETAIN - BETAI 

ERR1* CABS (TERROR) /CABS (BE TAIN) *100. 

ANG2 = AT AN ( A I MAG ( BFT A IN ) /RE AL (BET A IN) ) 

ERR2* ABS (ANG2 - ANG 1 ) 

ANG1 = ANG2 

IF (ID.GE.IDCR.AND.ERRl.GE.50.) 60 TO 540 
IF (ID.GE.IDMAX) GO TO 5*0 

IF (ID. EQ.l. OR. RK.LF. 0.0001) WRITE (6*665) ID.RETA IN.FRR1 *EPR? 
IF (ID.LT.5.0R.ERR1.GT.ERR0R.0R.ERR2.GT .ERROR) GO TO 60 
GO TO 545 
540 CONTINUE 

WRITE (6*365) 

545 CONTINUE 

WRITE (6*666) ID.BFTA1N*ERP1 »ERR2 
GO TO 550 

CHECK FOR CONVERGENCE ON OMEGA 

4000 CONTINUE 

CERROR s OMEGAN - OMEGA 
ANG2= PI/2. 

ERR 1 = 0.0 

IF (CABS(OMEGAN) .EO.O.O) GO TO 43 
ERRls CABS (CERROR) /CABS (OMEGAN) *100. 

IF (REAL (OMEGAN) .EO.O.O) GO TO 43 
ANG2= ATAN( A IMAG< OMEGAN) /REAL (OMEGAN) ) 

43 CONTINUE 

ERR2= ABS (ANG2-ANG1 ) 

ANG1 = ANG2 

IF ( ID.GE • I OCR. AND. ERR 1 *GE .50 • ) GO TO 547 
IF (ID.GE.IDMAX) GO TO 547 

IT (ID. EQ.l. OH. RK.LE. 0.0001 ) WRITE (6.666) ID. OMEGAN. ERR1 .ERR? 
IF (ID.LT.5.0R.ERR1.GT .ERROR. OR. ERR2.GT .ERROR) GO TO 60 
GO TO 546 
547 CONTINUE 

WRITE (6.365) 

546 CONTINUE 

WRITE (6*t>66) 10. OMEGAN. ERR1 .ERR2 


C 

C 

C 


C 

C 

C 


550 CONTINUE 

IF (ID.GE.IDMAX) WRITE(6.390> 

PRINT OUT FINAL MU MATRIX 

WRITE (6.2) 

WRITE (6.110) 

DO 106 Nisi, NTS 
Ns N1 - 1 

WRITE <6.:04> N. (MU(N1.U .Lsl.LTS) 
106 CONTINUE 
WRITE (6.2) 


CALCULATE AND PRINT EQUIVALENT BETAL IF CALCULATED FROM GEOMETRY 

IF (OPTION. GE. 3. AND. OPTION. LE. 5) WRITE (6,501) RO.AK1 
IF (OPTION. GE. 3. AND. OPTION. LE. 5) BETAL= BL INER/WA/PRES 
IF (OPTION. GE. 3. AND. OPTION. LE. 5) WRITE (6.370) BETAL 


**» >• 



non 


-48- 


C 

c 

c 


11 


19 


18 


CALCULATE AND PRINT EITHER N, TAU OR HETAI WHICHEVER IS APPROPRIATE 

IF (OPTION. NE.l. AND. OPTION. NE. 3) wRITF. (6.380) RETAIN 
IF (OPTION. EQ. 2. OR, OPTION. EQ. 4) GO TO 710 
WR = REAL (OMEGAN) 

WI = AIMAG(OMFGAN) 

BIR= REAL (BETAIN; 

611= AIHAG(BETAIN) 

INO= -((BIR - MACH/GAMMA) **2*HI 1**2) / (2. *MACH* (BIR-MACH/GAMMA) ) 

00 11 N=1 *3 

RN= FLOAT (N) - 1 . 

TAUO= ABS(ASIN(-8II/HACH/IN0) ♦ RN*PI)/WR 

BETAIN= MACH* (1. /GAMMA - IN0*(1. - CEXP (-CI*OMf GAN*TAUO) ) ) 

ERR1= ABS < (BIR - R C AL (BETAIN) ) /BIB) 

ERR2= ABS ( (B 1 1 - AIMAG(BETAIN) )/HlI) 

IF (ERR1. LE. 0.30. AND. ERR2.LE. 0.30) GO TO * 9 
RN= -PN 

TAU0 = A8S(ASIN(-BII/MACH/IN0) ♦ RN*PI)/WR 

BETAIN= MACH* (1. /Gamma - IN0*(1. - CEXP (-CI*OMEGAN*TAUO) ) ) 

ERR1= ABS ( (BIR - RFAL (BETA IN) ) /BIR) 

ERR2= ABS ( (Bl I - AIMaG (BETAIN) > /BI I ) 

IF (ERR1 »LE. 0.30. AND. ERW2.LF .0.30) GO TO 19 

CONTINUE 

WRITE (6*3) 

GO TO 710 
CONTINUE 

WRITE (6.386) RN. INO.TAUO.BFTAIN 
IF ( Wl.EU.0.0) GO TO 710 
RN= -3.0 
CONTINUE 
RN= RN ♦ 1.0 
IF (RN.GE.3.5) GO TO 17 
INI = INO 
T AU1 = TAUO 
DO IS 1=1.200 
TERM= bIR - MACH/GAMMA ♦ 

T AU= ABS (ATAN(-bI I/TERM) 


MACH* INI 
RN*P I ) /WR 


IN= TERM/(MACH*EXP(WI*TAU)*COS(WR*TAU) ) 

IF ( ABS ( T AU -TAU1) .LT. 1 .E-5.AND. ABS ( IN - I N 1 ) .L T . 1 .E-5 . ANO. I . OF . S ) 
1 GO TO 16 

IF (I. GE. 5. AND.IN.LT. 0.0) GO TO 16 
IF (I.GE. 5. AND.TAU.LT. 0.0) GO TO 16 
INI = IN 
T AU1 = TAU 

15 CONTINUE 

16 CONTINUE 

IF ( IN. LT. 0.00. OR. TAU. L r. 0.00) r,o TO 1H 

BETAIN= MACH* ( 1 ./GAMMA - IN*(1. - CEXP ( -C I *OMEGAN*TAU) ) ) 

ERR1 = ABS ( (BIR - RFAL (BETAIN) ) /bIR) 

ERR2= ABS ( (Bl I - AIMAG (BE TAIN) > /HI I ) 

IF (ERR1.GE.0.15.0R.FRR2.GE.0.15) GO TO 18 
WRITE (6.385) RN. IN. TAU.bETAlN 
GO TO 710 

17 CONTINUE 

SWT* SIN(WR«TAUO) 

CWT= COS(WR*TAUO) 

EWT= EXP <WI*TAU0) 

EC 1 = EWT*CWT - 1. 

WTW= WR/TAN(WR*TAUO) ♦ Wl 

TAU1 = (B IR/MACH/EC 1 ♦ BII/MACH/EWT/SWT - 1 . /GAMMA/E C 1 )/ (INO* 

1 (EWT* ( WI*CWT - WR*SWT ) /EC 1 - WTW) ) 

INI = -BII/MACH/EWT/SWT - INO - IN0*TAU1*WTW 
TAU= TAUO ♦ TAU1 
IN= INO ♦ INI 

BETAIN= MACH* ( 1 ./GAMMA - IN*(1. - CEXP ( -C I *OMEGAN*TAU ) ) > 

WRITE (6.5) IN. TAU. BETAIN 

CALCULATE RUNNING TIME FOR THIS SET OF DATA 


710 CONTINUE 

CALL SECOND 


(TENO) 


i 
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LINEARIZEO SOLUTION I' 
.F9.6.5H AND .F9.6./. 


TEX = TEND - TBGN 

WRITE (6.80) TBGN.TEND.TEX 

RETURN FOR NEXT SET OF OATA 

GO TO 570 
6000 CONTINUE 
STOP 

FORMAT STATEMENTS 

1 FORMAT ( 1H1 » /// ) 

2 FORMAT (///) 

3 FORMAT <48N NO AND TAU0 WOULD NOT CONVERGE FOP THIS PROBLEM,//) 

4 FORMAT (66H N AND TAU DID NOT CONVERGE SO A LINEARIZATION SOLUTI 
ION IS GIVEN,//) 

5 FORMAT (6IH N ANO TAU DID NOT CONVERGE SO A 

1 GIVEN. /.5X.29HTHE CALCULATED N AND TAU ARt 
25X.18HGIVING A BETAI OF ,F9.6.3X,F9.6,///> 

80 FORMAT (7H TBG=,F8.3,6K TEND=,F8.3,5H TEX=,F8.3,///) 

100 FORMAT <4F20.10»/,8F10.8./?7I10> 

101 FORMAT (8F10,8 ) 

104 FORMAT ( 2X. 13* 3X.8G 13,6. /.(8X.8G 13.6,/) ) 

105 FORMAT (49H THE FIRST GUESS FOR THE MU MATRIX IS AS F0LL0WS.//5X, 
11HN.16X,3HL=1.24X.3HL=2»24X»3HL=3./) 

110 FORMAT (35H THE FINAL HU MATRIX IS AS FOLLOWS. //5X . 1HN. 16X .3HL-1 , 
124X.3HL=2,24X,3HL=3,//) 

200 FOF.MAT ( 85H THE FOLLOWING CALCULATIONS ARE MADE FOR A COMBUSTOR WIT 
1H THE FOLLOWING CONFIGURATION.///) 

210 FORMAT (20X.23HTKE LENGTH TO RADIUS = .F9.6.24X.20HMEAN FLOW MACH N 
10. s,F9.6//.20X,7HLHAT = * 1 1 „5X, 7HMHAT = .I1.5X.7HNHAT =■• .I1.22X. 
29H0PTI0N = .I2.//.20X. 8HGAMMA = .F9„6,///> 

220 FORMAT (40H THE ABSORBER HAS THE FOLLOWING GEOMETRY .///) 

230 FORMAT (20X.20HTHE LINER STARTS AT .F9.6.13H AND ENDS AT .FQ.6.SX, 
117HAPERTURE WIDTH = ,F9.6,//) 

235 FORMAT (20X.8H8ETAL = .F9.6.3X.F9.6,///) 

240 FORMAT (20X.23HTHE BACKING DISTANCE = .F9.6.24X.L5HTHE APERTURE ME 
1 AN FLOW s ,F9.6,//,20X» 18H APERTURE LENGTH = ,F9.6v28X, 19HEFFECTI VE 

2 LENGTH = .F9.6.//.20X, 10HEPSILON = .F9.6.36X, 14HCAVITY WIDTH =♦ 

3 F9.6.//.20X, 16HCAVITY DENSITY = .F9.6.30X, 18HAPERTURE DFNSI TY = 

4 .F9.6.//.20X.20HCAVITY SOUND SPEED = .F9r6,///> 

250 FORMAT (39H THE INJECTOR RESPONSE IS DESCRIBED BY ,///) 

260 FORMAT (20X.8HBE7AI = ,F9.6»3X,F9.6,///'> 

270 FORMAT (20X.40HFOR THE CROCCO SENSITIVE TIME LA6 THEORY. 16X.4HN - 

1. F9.6.4X.6HTAU = .F9.6,// 

2 .20X.16H INITIAL BETAI = ,F9.6. 3X,F9.6»///> 

280 FORMAT (20X.80HFOR THIS OPTION FREQUENCY IS INPUT AND HFTAI IS Cu. 
1CULATED. INPUT FREQUENCY IS ,F9.6,3X.F9.6.//.20X,6HWWO = .F9.b,j< 

2. F9.6,/) 

285 FORMAT (20X.27HTHE FIRST GUESS OF BETAI IS.F9.6.3X.F9.6.///) 

290 FORMAT (37H THE NOZZLE RESPONSE IS DESCRIBED BY ,///) 

3C0 FORMAT (20X.43HTHE FUNDAMENTAL FREQUENCY FOR THIS MODF IS .F9.6.3X 
1.3H0.0,/) 

305 FORMAT (20X.27HTHE FIRST GUESS OF OMEGA IS.F9.6.3X.F9.6,///) 

310 FORMAT (20X.8HBETAN = ,F9.6,3X,F9.6,///) 

320 FORMAT (43H MISCELLANEOUS INFORMATION FOR THIS PROBLEM. ///.20X.26H 
1 NUMBER OF TERMS RADIAL .I2.4X. 1 3HLONGI TUDINAL .I2v9X,BHIDMAX = 

2.I4.//.20X.8HERROR * .F9.6,/) 

350 FORMAT (78H ITER REAL OMEGA IMAG OMEGA ERRORS 

1 MODULUS ANGLE.//) 

360 FORMAT (78H ITER REAL BETAI IMAG BETAI ERRORS 

1 MODULUS ANGLE.//) 

365 FORMAT ( IX ,5lH***THlS PROBLEM DOES NOT APPEAR TO BE CONVERGING***) 
370 FORMAT (27H THE EQUIVALENT BETAL IS .F9.6.3X.F9.6,///) 

380 FORMAT (21H THE FINAL BETAI IS.F9.6.3X.F9.6,///) 

385 FORMAT (28H N AND TAU CONVERGED FOR PN=.F9,6./» 

1 5X.29HTHE CALCULATED N AND TAU ARE .F9.6.5H AND .F9.6.16H 
2RESPECTIVELY./.5X.18HGIVING A BETAI OF ,F9.6,3X,F9.6.///) 

386 FORMAT >30H NO AND TAUO CONVERGED FOR RN=,F9.6,/, 

(THE CALCULATED NO AND TAUO ARE .F9.6.5H AND .F9.6.16H 

A BETAI OF .F9.6.3X.F9.6,///) 


1 5X.31H1 _ _ . 

2RESPF.CTIVELY./.5X, 18HGIVING 


CftGtflAL PAGE is 
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190 Fi.HMAT (/.44H tt «**» THIS PROBLEM HAS NOT CONVERGED # 0 ***,///) 

501 FORMAT <26H THE SLOT IMPFDANCES ARE ♦ 

1 10X.3HRO=.G13.6.5X.3HKI=,G13.6. 

2 ///) 

(S 66 FORMAT <2X.I3.bX,6l3.6.2X,Gn.6.15X.F10.6.5X.F10.6) 

900 FORMAT (SX.37HY0UR VALUE FOR MACH IS EXTREMELY HIGH . G1 0 . 4 ♦ // ) 

902 FORMAT (5X » 39HY0UR VALUE FOR LENGTH IS EXTREMELY HIGH. G] 0.4.//) 

904 FORMAT (5X.42HY0UR PRODUCT OF WA^BETAL IS EXTREMELY HIGH . ?G 1 0 . 4 , // ) 
906 FORMAT (SX.38HY0UR VALUE FOR 6FTAN IS EXTREMELY HIGH.2G1 0 .4 .// ) 

QOS FORMAT (SX.35HYOUR VALUE FOP BR IS EXTREMELY HIGH .G1 0 . 4 . // ) 

Q 1 0 FORMAT (SX.36HYOUR VALUE FOP AMF IS EXTREMELY HIGH . G1 0 .4 . // ) 

912 FORMAT (SX.35HY0UR VALUE FOP IN IS EXTREMELY HI6H.G10.4.//) 

Q 1 4 FORMAT (5X.36HY0UR VALUE FOP TAU IS EXTREMELY Hi GH .G1 0 . 4 . // ) 

916 FORMAT (SX.55HTHE MAGNITUDE OF YOUR VALUE FOR BETA I IS EXTREMELY H 
1 IGH.2G10.4.//) 

920 FORMAT ‘5X.37HY0UR VALUE FOP GAMMA IS EXTREMELY LOW .G1 0 .4 , // ) 

922 FORMAT (5X.38HY0UR VALUE FOR LFNGTH IS EXTREMELY LOW.G 1 0 . 4 * // ) 

924 FORMAT (5X.34HY0UR VALUE FOP BP IS FXTREMELY LOW.G10.4,//) 

926 FORMAT (SX.35HYOUR VALUE FOR AMF IS EXTREMELY LOW .61 0 . 4 . // ) 

930 FORMAT (5X.109HWITH THIS APERTURE THE INTEGRATION OVFR THE LlNF'R S 
lHOULO BE CARRIED OUT SEE COMMENT CARDS ABOVE 500 CONT INUF .G1 0 .4 , 
2//) 

940 FORMAT (5X.89HY0U HAVE NOT KEPT ENOUGH TERMS IN THE LONGITUOIMAl D 
1 IRECT I ON FOR THE MODE YOU HAVE CHOSEN. /.6HNTS = * I2.5X . 7 hNhAT = . 
212.//) 

942 FORMAT (5X.84HY0U HAVE NOT KEPT ENOUGH TERMS IN THE RADIAL DIRECTI 
ION FOR THE MODE YOU HAVE CHOSEN ./ .bX . 6HLTS = . I2.5X.7HLHAT = .12. 
2//) 

950 FORMAT (5X.61HTHIS FREQUENCY IS EXTREMELY HIGH FOR THE MODE YOU HA 
1VE CHOSEN.) 

960 FORMAT (5X.98HONE OF THE VALUES LISTED BELOW IS NOT PHYSICALLY Mp A 
1NINGFUL THIS SET OF DATA WILL NOT BE EXECUTED ./ *5X . 7HMACH = .G10. 
24.5X.8HGAMMA s .G10.4.5X.9HLENGTH s .G10.4.5X.BHBETAL = .2G10.4,/, 
35X.8HBETAN s . 2G1 0 .4 »5X .5HX A = .G10.4.5X.5HXB s .G10 .4 .5X .5HIN r , 
4G10.4.5X.6HTAU = .G10.4,////) 

980 FORMAT (SX.70HTHE ACOUSTIC MODE YOU SELECTED IS TOO HIGH FOR THIS 
1 PROGRAM TO HANDLE ./ .SX . 7HLHAT - ♦ 1 2 ,5X . 7HMHAT = • 1 2 »5X . 7HNHAT = « 

2 12) 

END 
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FUNCTION BESRT (M.L ) 

C* THIS FUNCTION SUBROUTINE IS A TABLE OF THE ROOTS OF THE RESSFL 

C* FUNCTION OF THE FIRST KINO. - 

DIMENSION A ( 4 * 1 0 ) 

DIMENSION B ( 1 f 1 0 ) 

IF (M.EG.O) GO TO 10 


A ( 1 

.1) 

= 

1.84118378 

A(1 

*2) 

= 

5.33144277 

A ( 1 

.3) 

S 

8.53631637 

Ad 

.4) 

= 

11.7060049 

A(1 

♦ 5) 

= 

14.86358863 

Ad 

.6) 


18.01552786 

A(1 

♦ 7) 

= 

21.16436986 

A(1 

.8) 

s 

24.31132666 

A(1 

♦ 9) 

s 

27.45705057 

A(l. 

10) 

= 

30.60192297 

A ( 2 

.1) 

= 

3.05423693 

A ( 2 

♦ 2) 

s 

6.70613319 

A ( 2 

♦ 3) 

s 

9.96946782 

A ( 2 

.4) 

s 

13.17037086 

A (2 

.5) 

s 

16.34752232 

A ( 2 

. 6 ) 

£ 

19.51291278 

A (2 

.7) 

£ 

22.67158177 

A ( 2 

.8) 

£ 

25.82603714 

A (2 

.9) 

£ 

28.97767277 

A (2. 

10) 

= 

32.12732702 

A ( 3 

.1) 

s 

4.20118894 

A ( 3 

.2) 

= 

8.01523660 

A ( 3 

*3) 

s 

11.34592431 

A ( 3 

.4) 

s 

14.505848 19 

A (3 

.5) 

s 

1 7 . 788747’ '7 

A<3 

.6) 

s 

20 • 972476 >4 

A ( 3 

*7) 

s 

24.14409743 

A (3 

.8) 

s 

27.31005793 

A ( 3 

.9) 

s 

30.47026881 

A (3. 

10) 

s 

33.62694918 

A (4 

*1) 

s 

5.31755313 

A (4 

*2) 

X 

9.28239629 

A (4 

» 3 ) 

s 

12.68190044 

A (4 

.4) 

s 

15.96410704 

A (4 

.5) 

s 

19.19602880 

A (4 

.6) 

r 

22.40103227 

A (4 

.7) 


25.58975968 

A (4 

.8) 

s 

28.76703622 

A (4 

♦ 9) 

s 

31.93853934 

A (4. 

10) 

= 

35.10391668 


BESRT=A(M,L) 

RETURN 
10 CONTINUE 

B ( 1 . 1 ) = 0.00000000 
B ( 1 ♦ 2 ) = 3.831 70597 
B ( 1 .3) = 7.01558667 
B ( 1 .4) = 10.17346814 
8(1 »5) = 13,32369194 
B ( 1 .6) = 16.47063005 
B ( 1 .7) = 19.61585851 
0 ( 1 .8) = 22.76008438 
0 C 1 .9) = 25.90367209 
B < 1 • 1 0 ) = 29.04602853 
BESPTs B(1,L) 

RETURN 

END 


i 
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C* THIS FUNCTION SUBROUTINE IS A TABLE OF THE VALUES OF THE BESSEL 
C* FUNCTION OF THE FIRST KIND. 

DIMENSION A ( 4 * 1 0 ) 

DIMENSION C(l*10) 

IF (M.EQ.O) GO TO 10 


A ( 1 

*n 

Z 

0.58186522 

A < 1 

.2) 

- 

-0.34612620 

A ( 1 

♦ 3) 

z 

0.27329994 

A ( 1 

.4) 

r 

-0.23330442 

A ( 1 

.5) 

z 

0.20701265 

A ( 1 

♦ 6) 

Z 

-0.18801749 

A ( 1 

♦ 7) 

Z 

0.17345905 

A ( 1 

.8) 

z 

-0.16183821 

A ( 1 

.9) 

z 

0.15228207 

A U t 

10) 

z 

-0.14424290 

A (2 

.1) 

z 

0.48649868 

M2 

.2) 

z 

-0.31353045 

A ( 2 

.3) 

z 

0.25474416 

M2 

♦ 4) 

z 

-0.22088158 

A (2 

.5) 

z 

0.19793743 

A (2 

.6) 

z 

-0.18101000 

A ( 2 

.7) 

z 

0.16783553 

A (2 

.8) 

r 

-0.15719517 

A 1? 

.9) 

z 

0.14836378 

A (2< 

10) 

z 

-0.14087833 

A (3 

*1 ) 

z 

0.4343942763 

A (3 

.2) 

z 

-0.2911584413 

A (3 

.3) 

z 

0.240738175 

A <3 

*4) 

z 

-0.210965204 

A (3 

.5) 

z 

0.190419022 

A (3 

.6) 

z 

-0.175048405 

A (3 

.7) 

z 

0.162954965 

A (3 

.8) 

z 

-0.153102409 

A C 3 

.9) 

z 

0.144866574 

A ( 3 1 

.10) 

z 

-0.137844513 

A (4 

.1) 

% 

0.3996514545 

A ( 4 

.2) 

z 

-0.2743809949 

A (4 

.3) 

z 

0.229590468 

A (4 

• 4) 

z 

-0.202763849 

A (4 

.5) 

z 

0.184029896 

A ( 4 

» 6 ) 

z 

-0.169878516 

A (4 

.7) 

z 

0.158655372 

A (4 

.8) 

z 

-0.149451156 

A (4 

.9) 

z 

0.141714307 

A (4 

.10) 

z 

-0.135086328 


RESVL=A(M.L) 

RETURN 
10 CONTINUE 

C(1 »D = 1.00000000 

C!1 .2) = -0.4027588095 
C(1 .3) = 0.301128303 

C(i .4) = -0.249704877 
C(1 .5) = 0.218359407 

C(1 .6) = -0.19645371 
C ! 1 .7) = 0.180063375 

C ( 1 .8) = -0.16 T 184600 
C < 1 *9) = 0.1 56724985 

C(l.lO) = -0.148011108 
RESVL= C<1.L> 

RETURN 

FND 
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42 


44 

45 
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SUBROUTINE LINER (XA*XB«PIL»L*N*MHAT»MUX»CIOM*MACH*BLlNER*LTS» 

*THIS C SUBROUTINE IS USED TO CARRY OUT THE INTEGRATION OVER THE 
LINER 

REAL NMNP* NPNP* MACH 

COMPLEX MUX * C I OM * BL I NER * SUM 1 *SUM2*CZERO 
DIMENSION MUX (50*10) 

BLINER= CZERO 
DO 47 LP=1.LTS 
SUM1= CZERO 
SUM2= CZERO 
DO 45 N1=1*NTS 
NP= N1 - 1 
RN= FLOAT (NP) 

RNP IL= RN*PIL 
IF (N.EQ.NP) GO TO 42 
NMNP= FLOAT (NP - N)*PIL 

SUMI= SUM1 T 1 N MUX (Nl!LPf<* (SIN (NMNP*XB> /2./NMNP ♦ SIN(NPNP*XB> /?./ 
1 NPNP - SIN(NMNP*XA)/2./NMNP - S IN ( NPNP* X A ) /2. /NPNP ) 

SUM2= SUM2 ♦ MUX (N1 » LP) *RNPIL* (COS (NMNP*XB) /2./NMNP ♦ COS(NPNP"» 

1 XB)/2./NPNP - COS (NMNP*XA)/2./NMNP - COS ( NPNP*X A ) /2 . /NPNP ) 

GO TO 45 
CONTINUE 

SUMl= P SUMi°i MUxInI^LP) ^ ( ( XB - XAJ/2. ♦ (SIN (2.*RNPIL*XR) - 

^uSiS f i*S5 lij CR W PIL-XB» **2 - SIN(RNPIL*XA>* 
1«2)/(2.*RNPIL> 

GO TO 45 
CONTINUE 

SUM1= MUX(N1«LP)«(XB - XA) 

RLINER= E BLINER - BFSVL(MHAT.LP>*<CIOM*SUMI ♦ MACH*SUM2> 

CONTINUE 
RETURN 
END 





